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Stability of FFLO states in optical lattices with bilayer structure 
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We investigate the stability of the superfiuid state in a bilayer fermionic optical lattice system with a 
confining potential, using the Bogoliubov de-Gennes equations. It is clarified that in the imbalanced case, 
the introduction of the interlayer hopping stabilizes the radial Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) 
state, while makes the angular FFLO state unstable. We also discuss the system size dependence of the 
superfiuid ground state. It is clarified that in a certain ring region the A-FFLO state is indeed realized in 
a large system. 
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1. Introduction 

Recently ultracold atomic gases have attracted much 
interest since the successful realization of Bose-Einstein 
condensation (BEC) in a bosonic ^''Rb system. Among 
them, an ultracold gas system in a periodic potential, so- 
called, an optical lattice system, has been providing 
an ideal stage for experimental and theoretical studies of 
fundamental problems in condensed matter physics. Due 
to its high controllability in interaction strength, par- 
ticle number, and other parameters, many remarkable 
phenomena have been observed such as a phase transi- 
tion between a Mott insulator and a superfiuid in bosonic 
systems,''^ and a crossover between the Bardeen-Coopcr- 
Schrieffer (BCS) state and the BEC state in fermionic 
systems. In addition, the superfiuid state in the 
spin-imbalanced fermionic systems has been realized in a 
fermionic ^Li system,^^' -'^^^ which stimulates further the- 
oretical investigations on the superfiuid state and its re- 
lated phenomena. 

An interesting question for the fermionic optical lat- 
tice system with imbalanced populations is how the 
symmetry-breaking state is realized at low tempera- 
tures. Various ground states have already been proposed 
to be more stable than the polarized superfiuid (PSF) 
state. ^"^'^^^ One of the probable candidates is the Fulde- 
Ferrell-Larkin-Ovchinnikov (FFLO) phase, ^^'^^^ where 
Cooper pairs are formed with a finite total momentum. 
This phase has been observed in the high field region in 
CeCoIus^^"-'^^^ and has theoretically been discussed in the 
compound^""^^' as well as cold atoms with imbalanced 
populations. ^''"•^"^ In the two dimensional optical lattice 
with a confining potential, it has been pointed out that 
two kinds of the FFLO states are realized such as the 
radial-FFLO (R-FFLO)3i'32) ^nd the angular-FFLO (A- 
FFLO) states. In the former state, the superfiuid or- 
der parameter changes its sign along the radial direction. 
In the latter, the order parameter oscillates along the an- 
gular direction, and the C^^ symmetry as well as U{\) 
symmetry are broken. On the other hand, such FFLO 
states with a three-dimensional structure have not been 
studied so well although the interlayer coupling should 



be important for realistic optical lattice systems. '^^^ 

To make this point clear, we systematically study 
the two-dimensional optical lattice system with a bi- 
layer structure as a simple model, which is schemati- 
cally shown in Fig. 1. We discuss how the interlayer hop- 




Fig. 1. (Color online) Optical lattice with a bilayer structure 

ping between two layers affects ground state properties 
by means of the Bogoliubov de-Gennes (BdG) equation. 
We then clarify how the R-FFLO and A-FFLO states are 
realized in the optical lattice with a confining potential. 
The system size dependence of the ground state is also 
discussed. 

The paper is organized as follows. In §2, we introduce 
the model Hamiltonian for the two-component fermions 
in the optical lattice and briefly summarize our theoret- 
ical approach. In §3, we discuss the stability of the R- 
and A-FFLO states in the bilayer optical lattice system. 
Scaling behavior in ground state properties is addressed 
in §4. A summary is given in the final section. 

2. Model and Method 

Wc study two-component ultracold fermions trapped 
in a two-dimensional bilayer optical lattice, which should 
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be described by the following Hubbard Hamiltonian, 



H = -t ^ [c\^^Cja„ + h.c. 

{ij)cro 
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where (ij) indicates the nearest neighbors in each layer 
with L X L sites. Ciaa {^laa) annihilates (creates) a 
fermion at the ith site of a{~ l,2)th layer with spin 
o" (=t , 4-), and h,aa = cj^aCiaa. t{t') is the intralayer 
(inter layer) hopping and U{< 0) is the attractive inter- 
action. The total number of particles N{=^ + N{) and 
the imbalanced population P[= {N^ - N^)/In^ + iVj,)], 

where = J2ia{4aa^iaa), are tuned by the chemical 
potential /i and the magnetic field h although these quan- 
tities can be controlled directly in the experiments. The 
confining potential is defined by Via = Pifia/ clY , where 
/?(> 0) is the curvature of the harmonic potential, is 
the distance measured from the center of the ath layer, 
and a is the lattice constant. 

In the mean-field approach, the interaction term 
should be given as 



^la. I 1 
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where A^q = U {ciiaCi^a) is the local superfluid order 



parameter, and rii, 



,) is the local particle 



density. The BdG equations should be written in terms 
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where H"^^ = 
—t'Sij, and Aij 



(3) 
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-t(5(jj) + (-^ - ha + Uniaff)6ij, K,j = 
- l^iSij, where 5ij is the Kronecker delta. 



The eigcnfunctions (/)^"'' = * ('"q"-/ 1 '^54 ) indicate the Bo- 
goliubov quasiparticle amplitudes at the ith site on the 
ath layer. The mean fields are then determined by the 
self-consistent equations as 
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where /(e) = [exp(e/T) + is the Fermi distribution 
function and T is the temperature. 

In the iterative method, one sometimes reaches the 
metastable state with the higher energy than the ground 
state. Therefore, it is important to choose an appropriate 
initial state in this treatment. In this paper, to discuss 
the stability of the FFLO states, we use the symmetry 
breaking states with n-fold symmetry in space and solve 
the BdG equations iteratively. Then we determine the 
ground state by comparing the free energies of their con- 
verged solutions. 

For convenience, we define the characteristic length of 
the potential as c? = \ftf^a and the effective particle 
density as p = No? . We set t as a unit of energy 
and fix the interaction strength and the temperature as 
U/t ^ -4 and T/t = 0.001, for simphcity. The temper- 
ature is low enough, which enables us to discuss ground 
state properties in the bilayer optical lattice system. The 
obtained ground state always has a mirror symmetry in 
the case t' /t ^ 0. Therefore, in this paper, we will show 
the results for one of the layers to discuss the stability of 
the superfluid state. 

3. Stability of the superfluid state 

In the section, we focus on the optical lattice system 
{L = 30) with the confining potential /3 = 0.04 [d = 5a) 
to discuss how the introduction of the interlayer hopping 
affects the stability of the BCS, R-FFLO and A-FFLO 
states. Before starting discussions, we would like to men- 
tion the characteristic profiles for these states. When the 
magnetic field is small enough, the BCS state is real- 
ized, where the finite pair potential appears without the 
magnetization. In the R-FFLO or A-FFLO states, the 
magnetization is finite, and the sign changes appear in 
the pair potential, depending on the direction of its os- 
cillation. Namely, the sign change in the radial direction 
appears in the former state. On the other hand, it ap- 
pears in a certain ring region in the latter state. Note 
that the BCS and R-FFLO states have the same spatial 
symmetry, in contrast to the A-FFLO state. Therefore, 
the BCS state is adiabatically connected to the R-FFLO 
state through the crossover. 

We first study ground state properties of the dilute 
system with p ~ 0.89 [N ~ 70). In the case with 
t' /t = h/t = 0, the particles are smoothly distributed 
and the order parameter appears around the center of 
the system. The local particle density Ni = riiaa 
and local order parameter A^ for each layer are shown in 
Fig. 2. When the magnetic field is applied, three kinds 
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Fig. 2. Profiles of the particle density Ni and the order param- 
eter Ai in the dilute system with p ~ 0.89, t' /t = 0.0 and P = 0. 



2 



J. Phys. Soc. Jpn. 



DRAFT 



of ground states appear. In the case h < h'{^ 0.7t), 
the BCS state is reahzed. In the intermediate region 
h' < h < hci^^ lAt), the R-FFLO state is reahzed, where 
the magnetization appears in the system and the order 
parameter changes its sign along the radial direction. Be- 
yond the critical field /ic, the pair potential is no longer 
finite, and the ground state is paramagnetic. Here we fo- 
cus on the system with h/t ^ 0.6 to discuss the effect 
of the interlayer hopping. The cross-sections of the or- 
der parameter and the magnetization, which is defined 
by rrii ~ riia^ — riia]^, arc shown in Fig. 3. When the in- 



phase diagram of the dilute system, as shown in Fig. 4. 
The crossover between the BCS state and the R-FFLO 
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Fig. 3. (Color online) Profile of the order parameter (a) and 
the magnetization (b) as a function of r/d in the dilute system. 



terlayer hopping t' is small, the system belongs to the 
BCS state, where the order parameter appears around 
the center of the potential (r/d < 1.7). The introduction 
of the interlayer hopping reduces the pair potential, typ- 
ically around r/d ^ 1.5. At last, the sign change appears 
in the superfluid order parameter although it may not be 
visible in the case t'/t = 0.2 in Fig. 3 (a). In addition, 
the magnetization is induced in the vicinity of the sign 
change point, as shown in Fig. 3 (b). This suggests the 
realization of the R-FFLO state. Further increase in the 
interlayer hopping yields a clear sign change in the or- 
der parameter and increases the magnetization monoton- 
ically. Therefore, we can say that the interlayer hopping 
stabilizes the R-FFLO state. 

By performing similar calculations, wc obtain the 




Fig. 4. Phase diagram of the dilute system with p ~ O.i 



state is roughly estimated by the appearance of the mag- 
netization, which is shown as a dashed line. The inter- 
layer hopping t' little affects the phase boundary be- 
tween the R-FFLO and normal states, in contrast to the 
crossover boundary. 

Next, we study ground state properties of the dense 
system with p ~ 6.4 {N ~ 500). In the case, a doubly 
occupied (Fock) state is realized around the center of the 
system due to the attractive interaction and the trap po- 
tential, as shown in Fig. 5 (a). Therefore, when h/t = 0.0, 
the superfluid state is realized away from the center and 
the doughnutlikc structure appears in the pair potential, 
as shown in Fig. 5 (b) . When the spin imbalanced popu- 




Fig. 5. Profiles of the particle density Ni and the order param- 
eter Ai in the dense system with p ~ 6.4, t'/t = 0.0 and P = 0. 



lation is introduced, the R-FFLO and A-FFLO states are 
realized in the cases {h' < h < hd) and [hd < h < hc2), 
where h' = 0.6t, hd = l.lt and hc2 = lAt. We here fo- 
cus on the latter case with h/t = 1.3, where oscillation 
behavior with nine peaks in the order parameter appear 
in a certain ring region, as shown in Fig. 6 (a). The in- 
troduction of the interlayer hopping t' monotonically de- 
creases the magnitude of the order parameters, as shown 
in Fig. 6. Finally it vanishes and the phase transition oc- 
curs to the paramagnetic state. This instability should 
be explained by the following. In the A-FFLO state, the 
oscillation of the order parameter appears in a certain 
ring region, which implies that one dimensional structure 
plays an important role to be stabilized. In this point of 
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Fig. 6. (Color online) Order parameter in the dense system 
when t' /t = 0.0, 0.6 and 1.0 (from the loft to the right). 



magnetization can be induced around the regions with 
Aia/t ^ 0. Therefore, when the large magnetic field is 
applied, many peaks appear in the profile of the order 
parameters. This tendency is essentially the same as the 
single layer optical lattice. Finally, the order parame- 
ter vanishes at /i = hc2 , where the phase transition occurs 
to the normal state. 

By performing the similar calculations, we obtain the 
phase diagram shown in Fig. 8. It is found that the inter- 



view, the interlayer coupling can be regarded as the in- 
troduction of the two-dimensional structure. Therefore, 
the A-FFLO state becomes unstable against the inter- 
layer coupling. 

To discuss how the magnetic field affects the nature of 
the superfluid state, we also show the number of peaks 
in the angular direction M and the average of the or- 
der parameter, which is defined by A = jAij/A'^, in 
Fig. 7. When h/t = 0.0, the BCS state is realized with 
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Fig. 7. (Color online) Squares and circles represent the average 
and the number of peaks in the angular direction of the order 
parameter in the dense system with t'/t = 0.6. 



K/t ^ 1.2 and M = 0. Applying the magnetic field, 
A is little changed up to /i = h'{^ 0.3i) and begins to 
decrease beyond it. This means that the crossover oc- 
curs to the R-FFLO state at h = h', where the sign 
change in the order parameter reduces A. Further in- 
crease in the magnetic field monotonically decreases the 
quantity. When h = hd, oscillation behavior suddenly 
appears in the angular direction of the order parame- 
ter (M = 6), and the phase transition occurs to the 
A-FFLO state. The R-FFLO and A-FFLO states have 
different spatial structures in the magnetization and or- 
der parameter. For example, the sign change in the or- 
der parameter appears in the radial and angular direc- 
tion, which means that these states are not adiabatically 
connected to each other around h = hd- This is con- 
trast to the results for one-dimensional system, where the 
second-order phase transition occurs between the BCS 
and FFLO states. '^^^ When the system approaches the 
phase boundary {hc2 ~ 1.4i), A is rapidly decreased and 
Af is increased. This may originate from the fact that the 




Fig. 8. Phase diagram for the dense system with p ~ 6.4 



layer hopping stabilizes the R-FFLO state, which is the 
same as the dilute case. On the other hand, the A-FFLO 
state becomes unstable against the interlayer hopping. 

We have treated the bilayer system to discuss how 
the interlayer hopping affects the stability of the FFLO 
states. Although the treated model is simple to discuss 
ground state properties of the layered optical lattice sys- 
tem, we believe that the obtained results capture the 
essence of the interlayer hopping. Namely, the A-FFLO 
state with interesting spatial properties survives in the 
system with a small interlayer hopping. In the following 
section, we discuss how the FFLO states are realized in 
the large system. 

4. Size dependence of FFLO states 

In the section, we study the stability of the FFLO 
states in the larger system. It is naively expected that 
low energy properties in the system are scaled by the 
characteristic length of the harmonic potential. By con- 
trast, the characteristic length of the FFLO state should 
depend on the spin imbalance in the system, which is 
not directly related to the length d, as discussed in the 
previous section. Therefore, it is necessary to clarify how 
interesting low energy properties are changed by the sys- 
tem size. 

To discuss scaling behavior in the system carefully, we 
fix the confining potential at four edges of the system 
as Via ~ 18t, and keep the effective particle density and 
the spin imbalance as p = 6.4 and P = 0.15. The cross- 
sections of the particle density and the order parameter 
are obtained, as shown in Fig. 9. In the BCS case (not 
shown), the profiles are well scaled, and thereby the state 
is stable in the large system. We find in Fig. 9 (a) that 
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(a) L=40 



(b) L=60 



(c) L=80 



in the R-FFLO state, a main peak structure in the order 
parameter is scaled. However, a size dependence appears 
around the regions with a sign change {r/d ~ 1.3 and 
2.2). Since the dip structures tend to shrink on the in- 
crease of the system size, it may be difficult to observe 
them in a large system, as a signature for the realization 
of the R-FFLO state. As for the A-FFLO state, disor- 
der behavior in the profiles of the order parameter and 
magnetization appears in the region (1.5 < r/d < 2.0), 
which is reflected by the oscillation of these quantities in 
the angular direction, as shown in Fig. 9 (b). Neverthe- 
less, their envelopes are well scaled, where the magnitude 
has a maximum at rg/d = 1.75. The spatial dependence 
of the order parameter and magnetization is also shown 
in Fig. 10. It is clearly found that each peak in the mag- 
netization nii is located at the sign change point along a 
certain ring in the order parameter A j . These are consis- 
tent with the fact that the A-FFLO state is realized in 
the ring region. Increasing the system size, we find that 
the number of peaks in a certain ring region is increased. 
This means that the peak structure in the angular direc- 
tion is not scaled. To clarify the size dependence of the 
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Fig. 9. (Color online) Profiles of the order parameter A^, mag- 
netization rui, and particle density Ni for the R-FFLO (a) and 
A-FFLO (b) states. 



oscillations, we also calculate the wave number defined by 
fcA-FFLO = M/ro, where M is the number of the peaks 
in the order parameter. Fig. 11 shows that it approaches 
a certain value A:a-fflo ^ 0.7/a when the system size is 
increased. These imply that there indeed exists a region 
with the A-FFLO state, and it is connected to the FFLO 
state expected in the one-dimensional Fermi gases. The 
characteristic wave number does not depend on the sys- 
tem size, but on the imbalanced populations. Therefore, 
we can say that the A-FFLO state should be realized 




Fig. 10. Order parameter (upper panel) and magnetization 
(lower panel) in the system with t'/t = 0.6, p ~ 6.4, P ~ 0.15, 
L = 40, 60, and 80 (from the left to the right). 



in the thermodynamic limit (L 
hopping is small enough. 
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Fig. 11. Squares (circles) represent the wave number character- 
istic of the A-FFLO state for the ground (meta-stablc) state. The 
number indicates peaks in the pair potential. 



Before closing this section, we wish to comment on 
the possibility of the supersolid state. In the attractive 
Hubbard model on the bipartite lattice in two or higher 
dimensions, it is known that the density wave ground 
state and the supcrfluid state are degenerate at half 
filling. Therefore, in a certain narrow region with 
(ui) ^ 1, the supersolid state, where both the density 
wave state and superfluid state coexist, may be real- 
ized. It is necessary to carefully deal with particle cor- 
relations beyond the BdG mean-field approach, which is 
now under considerations. 

5. Summary 

We have studied the stability of the superfluid state 
in a bilayer fermionic optical lattice system with imbal- 
anced populations by means of the BdG equations. It has 
been clarified that the introduction of the hopping be- 
tween two layers stabilizes the radial FFLO state, while 
makes the angular FFLO state unstable. We have also 
discussed scaling behavior in the superfluid state. It has 
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been found that the wave number characteristic of the 
A-FFLO state approaches a certain value when the sys- 
tem size is increased. This imphes that the A-FFLO state 
should be realized in the thermodynamic limit. 
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